%% Stock Return Extrapolation, Option Prices, and Variance Risk Premium
% By Adem Atmaz
% This file generates Table 7: Effects for the variance risk premium: Model vs. evidence
% It also requires 1 function file: olshac.m
% Typically it takes less than 10 seconds to run

clear all;  clc;
format compact
format

% Parameter Values follow from  "Atmaz_Table_4_ParameterValues.m"
alpha=0.50;
theta=0.25;
mu=0.0212;
kappa=0.3567;
zeta=0.0044;
sigma=0.1625;
rho=0;
r=0.0056;

D0=100;
a=2.75;
vec_theta=[0, 0.25, 0.50, 0.75];
N_th=length(vec_theta);

% Initilization of vectors
vec_RV=zeros(1,N_th);
vec_VSR=zeros(1,N_th);
vec_VRP=zeros(1,N_th);

for i=1:N_th
    theta=vec_theta(i);    % Vary theta
    % Implied Quantities after fixing theta
    M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
    M_c=theta/(alpha*(1+theta));
    M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
    M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
    M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
    M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
    M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
    M_corr_sv=M_Cov_con/sqrt(M_Var_con);
    
    % VSR and VRP quantities
    M_zeta_v=zeta*M_Var_con;
    M_kappa_v=kappa+(sigma*M_Cov_con);
    M_eta_v=(theta/M_b)*(1-(alpha*M_c))*M_Var_con;
    M_sigma_v=sigma*sqrt(M_Var_con);
    
    % long-run mean (steady-state) values of processes
    V0=zeta/kappa; % long run mean of  fundamental  variance
    v0=M_zeta_v/kappa;  % LRM of stock return variance v
    M_RN_LRM_vt=(M_zeta_v+(r*M_eta_v))/(M_kappa_v+(0.5*M_eta_v));
    X0=(r+(0.5*M_Var_con*V0))/(1+theta); % LRM of sentiment X
    A_Var_v=((M_sigma_v^2)/(2*kappa))*v0; % Var[v]
    temp1=(0.5*alpha*v0)/(1+theta);
    temp2=(((M_sigma_v^2)/(4*kappa))+(M_corr_sv*M_sigma_v))/((alpha*(1+theta))+kappa);
    A_Var_X=temp1*(1+temp2);
    A_Cov_vX=alpha*v0*temp2;
    
    tt1=M_kappa_v+alpha;   % temporary constant
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  #######################      instantenous  RV       ###################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    RV_A= M_zeta_v;
    RV_B=1-kappa;
    vec_RV(i)=(RV_A+(RV_B*v0))*(10000/12); %normalized to get monthly values
    
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  #######################   instantenous     VSR      ###################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    VSR_A=M_zeta_v;
    VSR_B=1-M_kappa_v;
    VSR_C=M_eta_v;
    vec_VSR(i)=(VSR_A+(VSR_B*v0)+(VSR_C*X0))*(10000/12); %normalized to get monthly values;
    
end

% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%  #######################     instantenous   VRP       ##################
%  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
vec_VRP=vec_RV-vec_VSR;

% REPORT MOMENTS
RV=round(vec_RV, 2)
VSR=round(vec_VSR,2)
VRP=round(vec_VRP,2)



disp('Press any key to continue to Predictive regression ...');
pause;

tau=1/12;  % variance horizon in years
h=3;       % predictive regression holding return horizon in months 1,3,12
N_sim=1000;                             % number of simulations: use 1000
S_end_year=18;                          % sample size 18 yeears or 62 years
int = 1/252;                            % time interval for each year in a path
time_grid = 0:int:S_end_year;           % Time vector from 0 to end year
N_time=length(time_grid);
N_rets=N_time-1;                        % Number of returns in each path
H=1/12;     % Data frequency in the regression e.g. 1/12 for monthly
N_pp=H/int; % number of data points for each data point in regression: 21 for month
% omitting the first burn out period
keep_final_years=S_end_year; % keep all of it
Obs=keep_final_years*12;

% 1 path vectors
vec_D = zeros(1,N_time); % Allocate output vector
vec_V = zeros(1,N_time); % Allocate output vector
vec_X = zeros(1,N_time); % Allocate output vector

% All paths monthly matrices
N_month=(12*S_end_year)+1;

SM_V = zeros(1,N_month); % Allocate output vector
SM_X = zeros(1,N_month); % Allocate output vector
SM_D = zeros(1,N_month); % Allocate output vector
SM_S = zeros(1,N_month); % Allocate output vector
SM_VSR = zeros(1,N_month); % Allocate output vector
SM_VRP = zeros(1,N_month); % Allocate output vector
SM_RV  = zeros(1,N_month); % Allocate output vector

% Univariate regression of VRP
vec_beta_VRP        =zeros(N_sim,N_th);
vec_R2_adj_VRP      =zeros(N_sim,N_th);
vec_SE_hac_VRP      =zeros(N_sim,N_th);
vec_tstat_hac_VRP  =zeros(N_sim,N_th);
vec_STD_U_VRP        =zeros(N_sim,N_th);
vec_STD_U_Ret        =zeros(N_sim,N_th);

% Univariate regression of RV
vec_beta_RV        =zeros(N_sim,N_th);
vec_R2_adj_RV      =zeros(N_sim,N_th);
vec_SE_hac_RV      =zeros(N_sim,N_th);
vec_tstat_hac_RV   =zeros(N_sim,N_th);
vec_STD_U2_RV        =zeros(N_sim,N_th);

% Joint regression of VRP and RV
Joint_vec_beta_VRP    =zeros(N_sim,N_th-1);
Joint_vec_beta_RV     =zeros(N_sim,N_th-1);
Joint_vec_SE_VRP_hac    =zeros(N_sim,N_th);
Joint_vec_SE_RV_hac     =zeros(N_sim,N_th);
Joint_vec_tstat_VRP_hac =zeros(N_sim,N_th);
Joint_vec_tstat_RV_hac  =zeros(N_sim,N_th);
Joint_vec_R2_adj        =zeros(N_sim,N_th-1);


% Begin simulation and regression
rng('default'); % fix the random number generator  seed  for replication
vec_Z1=randn(N_sim,N_time);
vec_Z2=randn(N_sim,N_time);
tic
for s = 1:N_sim
    s
    for th=1:N_th
        theta=vec_theta(th);
        % Stock price quantities
        M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
        M_c=theta/(alpha*(1+theta));
        M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
        M_a=a;
        M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
        M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
        M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
        M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
        
        % VSR and VRP quantities
        M_zeta_v=zeta*M_Var_con;
        M_kappa_v=kappa+(sigma*M_Cov_con);
        M_eta_v=(theta/M_b)*(1-(alpha*M_c))*M_Var_con;
        v0=M_zeta_v/kappa;
        M_RN_LRM_vt=(M_zeta_v+(r*M_eta_v))/(M_kappa_v+(0.5*M_eta_v));
        tt1=M_kappa_v+alpha;
        kappa1=(0.5*tt1)-(0.5*sqrt((tt1^2)-(4*alpha*(M_kappa_v+(0.5*M_eta_v)))));
        kappa2=(0.5*tt1)+(0.5*sqrt((tt1^2)-(4*alpha*(M_kappa_v+(0.5*M_eta_v)))));
        EE0=(1-exp(-kappa*tau))/(kappa*tau);
        EE1=(1-exp(-kappa1*tau))/(kappa1*tau);
        EE2=(1-exp(-kappa2*tau))/(kappa2*tau);
        
        % RV A and B
        RV_A=v0*(1-EE0);
        RV_B=EE0;
        % VSR A and B and C
        VSR_A=(M_RN_LRM_vt*(1-EE1-((kappa1/(kappa2-kappa1))*(EE1-EE2))))...
            +((M_zeta_v/(kappa2-kappa1))*(EE1-EE2));
        VSR_B=EE1-(((M_kappa_v-kappa1)/(kappa2-kappa1))*(EE1-EE2));
        VSR_C=(M_eta_v/(kappa2-kappa1))*(EE1-EE2);
        % VRP A and B and C
        VRP_A=(v0*(1-EE0))-VSR_A;
        VRP_B=EE0-VSR_B;
        VRP_C=-VSR_C;
        
        % Initalization of each path
        D0=100;
        V0=zeta/kappa; % long run mean of  fundamental  variance
        X0=(r+(0.5*M_Var_con*V0))/(1+theta); %  inital  X is equal to its long run mean
        vec_D(1)=D0;
        vec_V(1)=V0;
        vec_X(1)=X0;
        
        for j = 1:N_rets
            Z1=vec_Z1(s,j);
            Z2=vec_Z2(s,j);
            vec_V(j+1) = max(0, vec_V(j))+((zeta-(kappa*max(0, vec_V(j))))*int)...
                +(sigma*rho*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(sigma*sqrt(1-rho^2)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_X(j+1) = vec_X(j)+((r+(0.5*M_Var_con*vec_V(j))-(1+theta)*vec_X(j))*alpha*int)...
                +(alpha*M_sigma_1*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(alpha*M_sigma_2*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_D(j+1) = vec_D(j)+((mu*vec_D(j))*int)+(vec_D(j)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1);
            
        end
        
        % Get end of month values for each process from daily values
        SM_V= vec_V(1:N_pp:end);
        SM_X= vec_X(1:N_pp:end);
        SM_D= vec_D(1:N_pp:end);
        SM_S= SM_D.*exp(M_a+(M_b*SM_V)+(M_c*SM_X)); %stock price
        
        
        % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        %  #######################      instantenous   RV     ####################
        %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        RV_A= M_zeta_v;
        RV_B=1-kappa;
        
        
        % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        %  #######################   instantenous      VSR      ##################
        %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        VSR_A=M_zeta_v;
        VSR_B=1-M_kappa_v;
        VSR_C=M_eta_v;
            
        % simulated VSR VRP and RV
        SM_RV=(RV_A+(RV_B*M_Var_con.*SM_V))*(10000/12); %normalized to get monthly values in percentage squared
        SM_VSR=(VSR_A+(VSR_B*M_Var_con.*SM_V)+(VSR_C.*SM_X))*(10000/12);  %normalized to get monthly values in percentage squared
        SM_VRP=SM_RV-SM_VSR;       
        
        % Keep only the selected final years and omit the first years as burnout if needed
        SM_RV=SM_RV(end-Obs:end);
        SM_VRP=SM_VRP(end-Obs:end);
        SM_S=SM_S(end-Obs:end);
        SM_D=SM_D(end-Obs:end);        
        
        % Get h period ahead log returns
        SM_Ret_h=log(SM_S(1+h:end))-log(SM_S(1:end-h)); % log ret
        %SM_Ret_h=(SM_S(1+h:end)./SM_S(1:end-h))-1;      % returns
        SM_Ret_h_ann=(12/h)*SM_Ret_h*100;  %annualized in percentage terms
        
        
        % Univariate Regression of h period returns on VRP
        % ===============================================        
        reg_Y=SM_Ret_h_ann';            %   annualized   returns
        LT=length(reg_Y);
        reg_X=[ones(LT,1) SM_VRP(1:end-h)'];
        lag=round(0.75*(LT^(1/3)));    %  recommended hac lag choice        
        vec_STD_U_Ret(s,th)=std(reg_Y);  % get std for standardization later on
        vec_STD_U_VRP(s,th)=std(SM_VRP(1:end-h));
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        
        vec_beta_VRP(s,th)=beta(2);
        vec_tstat_hac_VRP(s,th)=t_NW(2);
        vec_R2_adj_VRP(s,th) = R2adj;
        
        
        % Univariate Regression of h period returns on RV
        % ===============================================        
        reg_X=[ones(LT,1) SM_RV(1:end-h)'];
        vec_STD_U2_RV(s,th)=std(SM_RV(1:end-h));
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        
        vec_beta_RV(s,th)=beta(2);
        vec_tstat_hac_RV(s,th)=t_NW(2);
        vec_R2_adj_RV(s,th) = R2adj;
        
        
        if theta > 0   % joint regression except for benchmark model theta=0
            
            % Joint Regression of h period returns on VRP and RV
            % ===============================================
            
            reg_X=[ones(LT,1) SM_VRP(1:end-h)' SM_RV(1:end-h)'];            
            [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
            
            Joint_vec_beta_VRP(s,th)=beta(2);
            Joint_vec_beta_RV(s,th)=beta(3);
            Joint_vec_tstat_VRP_hac(s,th)=t_NW(2);
            Joint_vec_tstat_RV_hac(s,th)=t_NW(3);
            Joint_vec_R2_adj(s,th) = R2adj;
            
        end  % end if for joint regression
        
    end   % end theta
    
end    % end sim
clc;
toc


%  REPORT   UNIVARIATE REGRESSION on VRP QUANTITIES
STD_beta_VRP=round(mean((vec_beta_VRP.*vec_STD_U_VRP)./vec_STD_U_Ret),2)
tstat_hac_VRP=round(mean(vec_tstat_hac_VRP),2)
 
%  REPORT  UNIVARIATE REGRESSION on RV  QUANTITIES
STD_beta_RV=round(mean((vec_beta_RV.*vec_STD_U2_RV)./vec_STD_U_Ret),2)
tstat_hac_RV=round(mean(vec_tstat_hac_RV),2)
 
%  REPORT   JOINT REGRESSION QUANTITIES
Joint_STD_beta_VRP=round(mean((Joint_vec_beta_VRP.*vec_STD_U_VRP)./vec_STD_U_Ret),2)
Joint_tstat_VRP_hac=round(mean(Joint_vec_tstat_VRP_hac),2)
Joint_STD_beta_RV=round(mean((Joint_vec_beta_RV.*vec_STD_U2_RV)./vec_STD_U_Ret),2)
Joint_tstat_RV_hac=round(mean(Joint_vec_tstat_RV_hac),2)
 

clear alpha D0 EE0 EE1 EE2 EEv j jj kappa kappa1 kappa2 m M_a M_b M_c;
clear M_corr_sv M_Cov_con M_DELTA M_eta_v M_kappa_v M_Ob_LRM_vt M_r M_RN_LRM_vt;
clear M_sigma_1 M_sigma_2 M_sigma_v  M_Var_con M_zeta_v mu n rho sigma;
clear theta tt1 v V0 Z1 Z2 zeta;
clear vec_Z1 vec_Z2;
 
%% THE END










